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We present an exact algorithm for finding all the ground states of the two-dimensional Edwards- 
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change with increasing system size and and with increasing antiferromagnetic bond ratio x. We find 
that that some system properties have very large and strongly non-Gaussian variations between 
realizations. 
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I. INTRODUCTION 



The Edwards- Anderson (EA) ±J spin glassb is a 
canonical example of a system with competing in- 
teractions that give rise to large numbers of low- 
energy states. Despite extensive investigation, the 



low-energy landscape 
versialr-ui-. hatb l _ l tb. 



sion 




still contro- 
two dimen- 



Although the two-dimensional ±J Edwards-Anderson 
model is much simpler than a three-dimensional case (it 
does— not have a phase transition at nonzero tempera- 
turet3 and individual ground states can be found kj-,a 
time that scales as a polynomial of the system sizce3), 
it still has many metastable states and a complex en- 
ergy landscape. At low temperatures spin glass relax- 
ation times become very long, complicating investiea- 
tions using standard Monte Carlo sampling techniquesEH, 
and also to varying extents more sauhiaticated sampling 



methods such_as, chaster algorithmslilrc-j and multicanon- 
ical methodsE 



jester algorithmsEJ't 

For ground state properties, exploiting optimization 
algorithms that find-exact, ground states has proven a 
powerful approachE3ll§E3H. However, these algorithms 
find a single ground state of a single realization, and one 
must sample appropriately from the ground states of each 
realizationEloE3 and also perform a reliable realization 
average to obtain correct results. 

In this paper we investigate the low-energy properties 
of the two-dimensional ±J E-A model by finding exactly 
all the ground states of each realization. We check that 
our enumeration is exhaustive by comparing the number 
of ground states that are found to exact results for the 
partitioa-function obtained using the method of Saul and 
KardaieH 

Though our algorithm is based on an existing 
polynomial-time algorithm that finds individual ground 
statesE3, it does not run in polynomial time. This is im- 
possible because the time just to enumerate the ground 
states grows exponentially with system size. Nonetheless, 
the number of ground states is vastly smaller than the 
number of spin configurations, and empirically we find 
that our algorithm runs in a time roughly linear in the 
number of ground states. Memory issues limit our cur- 
rent implementation to about 2 x 10 6 ground states. Be- 



cause there are huge variations in the number of ground 
states among realizations, the system sizes that we inves- 
tigate are rather small. Though the median number of 
ground states of a 10 x 10 system in which half the bonds 
are antiferromagnetic is 10 4 , at this system size 3 percent 
of the realizations have greater than 2x 10 6 ground states. 
Therefore, most of the data presented here is for systems 
of size 10 x 10 and less. 

The advantage of our method is that it produces qual- 
itatively new information because all the ground states 
are known explicitly and exactly, so that one can compute 
in detail the relationships between them. Moreover, these 
exact results can be used to validate sampling methods 
appropriate for larger systems. 

We use our method to investigate the paramagnetic- 
ferromagnetic phase transition that occurs as x, the 
fraction pf.. ajptiferromagnetic bonds in the system, is 
increasedll3E3ca. Quantitative analysis is complicated 
greatly by the fact that many quantities exhibit large, 
non-Gaussian variability between realizations. Nonethe- 
less, we are able to obtain new insight into the nature 
of the paramagnetic-ferromagnetic phase transition that 
occurs as the fraction of antiferromagnetic bonds is in- 
creased. 

The paper is organized as follows. Sec. II presents the 
algorithm, Sec. HI presents data on the distribution of 
presents results on the algorithm 



ground states, Sec. IV 
performance, and Sec. 



pre: 

[yj describes our investigation of 
the destruction of ferromagnetic order as the fraction of 
antiferromagnetic bonds is increased. The results are dis- 
cussed in Sec. VI. The appendix [A] gives a detailed pre- 



sentation of the algorithm. 



II. MODEL AND METHODS 



The Edwards- Anderson model 



(EA) modellil, in which nearest-neighbor Ising spins (ai — 
±1) on an L x L square lattice interact either via a ferro- 
magnetic or an antiferromagnetic coupling. The Hamil- 
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FIG. 1: Two ground states of a frustrated plaquette with 
four bonds. Ferromagnetic bonds are black lines, while anti- 
ferromagnetic bonds are dashed lines. Every configuration of 
spins produces at least one unsatisfied bond (denoted with a 
line perpendicular to the bond), and there are four minimum 
energy configurations. 



tonian is 



(1) 



where the sum (ij) is over all pairs of nearest-neighbor 
spins. Each bond is chosen to be cither +1 (fer- 
romagnetic) or —1 (antiferromagnetic) . We designate 
the fraction of antiferromagnetic bonds = is 

the Ising ferromagnet (no disorder), x = .5 (with equal 
numbers of ferromagnetic and antiferromagnetic bonds) 
is the maximally-frustrated spin glass (maximum disor- 
der), and x = 1 is the Ising antiferromagnet (no disor- 
der). Our systems range from x = .05 to x = .5 and have 
periodic boundary conditions. 



B. Calculating all the ground states of the EA 
model 



Our algorithm for finding all the ground states of the 
EA model first converts the problem of finding gr ound 
states into a graphical matching problem, as in Refs.E£lE2l. 
Next, all possible optimal matching solutions of this 
problem are found, and finally, these matchings are con- 
verted back into spin configurations. 



1. Conversion of energy minimization to a matching 
problem 



RefsHj'tJ show that the problem of finding a ground 
state for this spin glass model can be converted to a 
matching problem in graph theory. Here, we sketch out 
this conversion and discuss some subtleties that arise 
from our use of periodic boundary conditions. 

A ground state of a spin glass can be described not 
only in terms of spins and bonds, but-also as frustrated 
plaquettes and paths of broken bondstj. In a frustrated 
system, it is ipt possible for all bonds to be satisfied 
simultaneousljfij, which leads to a natural degeneracy of 
states. A simple example is shown in Figure [l]. 
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FIG. 2: Sample ground state of a spin glass, showing frus- 
trated plaquettes, unsatisfied bonds, and the corresponding 
spin and bond configuration. 




FIG. 3: Correspondence of spin glass ground state to solution 
of matching problem. The numbers along each edge of the 
matching solution indicate the weight of that edge. 



We denote plaquettes with an odd number of unsatis- 
fied bonds as frustrated, while satisfied plaquettes have 
an even number of unsatisfied bonds. Frustrated bonds 
form paths that connect frustrated plaquettes to each 
other. Because every frustrated plaquette has an odd 
number of frustrated bonds, it must be the endpoint of 
a path. Satisfied plaquettes either have no path through 
them or are midpoints in a path. This can be seen in 
Figure ||, where perpendicular lines have been added to 
frustrated bonds to show the paths. 

We identify the frustrated plaquettes as nodes of a 
graph, and the paths as edges with a weight equal to 
the number of broken bonds along the path. Ground 
states have the minimum number of frustrated bonds, so 
the problem of finding a spin glass ground state is also 
the problem of finding those edges that have the shortest 
total length. This problem arises in the context of graph 
theory and is. called the minimum weight perfect match- 
ing problem?]!. In solutions of this problem, each node is 
joined to one and only one other node, with the smallest 
possible total weight (which corresponds to the lowest 
possible energy). Figure ^ illustrates a sample ground 
state and its equivalent matching solution. 



2. Boundary conditions 

RefsE9'EH discuss the relation between spin glass 
ground states and solutions to a graphical matching prob- 
lem and prove a number of results for planar graphs 
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(a) Periodic (b) Antipcriodic 

NSEW NS, Periodic EW 




(c) Periodic NS, (d) Antiperiodic 

Antiperiodic EW NSEW 



FIG. 4: The four different boundary conditions for a sample 
lattice. Ferromagnetic bonds are solid lines and antiferromag- 
netic bonds are dashed lines. Bonds on the right and bottom 
sides wrap around the lattice to reconnect on the other side. 
Each boundary condition has the same set of frustrated pla- 
quettes, shown as filled-in circles. 



(e.g. free boundary conditions). In these references, the 
ground state problem is first converted to a, matching 
problem that can be solved in polynomial timcclrEa. This 
matching solution is then shown to correspond always to 
a spin glass ground state. 

For a periodic lattice, the transformation from spins 
and bonds to nodes and edges proceeds exactly as in the 
planar case, and the resulting matching problem can be 
solved in polynomial time. The issue that distinguishes 
this problem from the planar case is the conversion of 
the matching solution back to a ground state solution. 
The matching solution found will not always correspond 
to a ground state spin configuration for a given toroidal 
boundary condition. This complication arises because 
four lattices with four different boundary conditions will 
produce the same matching problem. These four bound- 
ary conditions are periodic on all sides, antiperiodic on 
the top and bottom, antiperiodic on the left and right, 
and antiperiodic on all sides. 

Boundary conditions can be changed from periodic to 
antiperiodic on an L x L system either by setting J™£ w = 
— JiL along the desired edge or by flipping the spins so 
that S , " 1 etu = —Sii along the desired edge. In this study, 
we flip the bonds. A sample lattice with the four different 
boundary conditions is shown in Figure 0. 



These four lattices have exactly the same frustrated 
plaquettes, so they produce the same matching problem. 
In this sense, the matching problem does not understand 
boundary conditions. When a matching solution is con- 
verted back into spins and bonds, it may not correspond 
to a ground state for a given boundary condition. 

We resolve this ambiguity by converting explicitly each 
matching into a spin configuration and checking the via- 
bility of each spin configuration for each boundary condi- 
tion. A ground state is only accepted for a given bound- 
ary condition if it has a consistent spin configuration. We 
find numerically that a matching solution always corre- 
sponds to a ground state solution of at least one boundary 
condition. 

This subsubsection has described the necessary pro- 
cedures for generating all the ground states in the case 
of periodic boundary conditions. The algorithm works 
perfectly for planar graphs without these procedures. 

3. Generating all optimal matchings 

Our algorithm for finding all the optimal matchings 
has three parts. The first part finds all edges that make 
up the optimal solutions. This part exploits the struc- 
ture of the edges, since the number of edges appearing 
in the ground states is a small subset of the total num- 
ber of edges. The second part takes this subset of edges 
and combines them to find all optimal matchings. The 
third part converts the optimal matchings into ground 
state configurations. The next section describes our al- 
gorithm briefly. More detail, including an example, is in 
the appendix |A|. r-ir-i 

The algorithm uses Edmonds' blossom algorithincilli3. 
which finds a single optimal solution to a matching prob- 
lem in polynomial time. W e use the Concorde implemen- 
tation of the algorithmtOcj. 

C. Finding all edges in all solutions 

The algorithm begins by making a list of nodes and 
possible edges. All frustrated plaquettes are found and 
designated as nodes. Pairs of nodes that are within a dis- 
tance of five are considered to have edges between them. 
This restriction controls the combinatorial explosion of 
possible edges, and optimal solutions involving weights 
larger than five are incredibly rare. Moreover, any errors 
introduced by this truncation are identified and elimi- 
nated at a later stage when the total number of ground 
states found for a realization is compared to an indepen- 
dent determination of the ground state degeneracy. 

To construct a list of edges that exist in at least one 
minimal weight matching, which we designate as viable 
edges, we start with an empty list. The blossom algo- 
rithm is run on the unmodified matching problem and 
a single matching solution is found. Each edge in this 
solution is added to the list of viable edges. 



To find more viable edges, the nodes are considered 
successively. For each node, a modified list of edges is 
created from the original list by deleting those known 
viable edges connecting to the current node. The blossom 
algorithm is run on this modified list to find an optimal 
solution for this new problem. If the solution has the 
correct path length (i.e. corresponds to a ground state), 
then the new viable edges that have been found are added 
to the list of viable edges. The process continues for this 
node. If the path length of the new solution is too large 
(i.e. does not correspond to a ground state), then we 
know that we have found all the viable edges associated 
with this node. The algorithm then proceeds to the next 
node in the list. By moving sequentially through the 
nodes, all viable edges are found. 



D. Determining optimal matchings 

The next part of the algorithm uses the list of viable 
edges to find all of the optimal matchings. It picks edges 
systematically from the list of viable edges until each 
node is connected by a given edge to one and only one 
other node. All possible combinations of viable edges 
in which each node is incident on exactly one edge are 
examined. 

Whenever there is this kind of combination of elements, 
there is a danger of a combinatorial explosion. In this 
case, the number of matchings (combinations of edges) 
is relatively controlled. Sec. IV discusses this issue in 
detail. 



E. Converting matchings to ground states 

All optimal solutions to the matching problem must 
be converted back into ground state spin configurations. 
This conversion is nontrivial because one matching so- 
lution can correspond to many different ground states, 
and the same ground state can be represented by differ- 
ent matchings. Simple examples of this phenomenon can 
be seen in Figure |5|. To resolve these complications, we 
keep every ground state we find in memory. Any pro- 
posed ground state is checked to see that it does not 
correspond to a ground state we have already found. 

The other important issue is the role of boundary 
conditions in this conversion from matchings to ground 
states. We determine the ground state or ground states 
from the matching by fixing the value of a single spin (in 
our case, we fix the upper left-hand spin as +1). Every 
other spin follows from this initial spin, because we know 
the specific bonds of the current boundary condition and 
their status as satisfied or unsatisfied. 
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FIG. 5: The relationship between matching solution is not 
one-to-one. As shown in the top diagram, two or more dif- 
ferent ground states can correspond to the same matching 
solution. In addition, as seen in the bottom diagram, a single 
ground state can correspond to multiple matching solutions. 
The filled circles are frustrated plaquettes, the thick dark lines 
are unsatisfied bonds, the open circles are nodes, and the thin 
dark lines are edges with lengths as shown. 



F. Partition function 

To check that our algorithm finds every ground state, 
we also generate the partition function of the realization 
at T = 0. This partition function gives the number and 
energy of the ground states of a given realization. We 
generate a partition function i.n polynomial time using 
Saul and Kardar's technique ,Edc3 which is a generaliza- 
tion of methods used for finding the partition fiinptinn 
at T = for the two-dimensional Ising model. ac3o - cZl 
For reasons of computational efficiency, we consider only 
LxL lattices where L is even. Because our methods yield 
ground states not only for lattices with regular periodic 
boundary conditions but also for those with antiperiodic 
boundary conditions, we generate four different partition 
functions for each possible lattice, corresponding to the 
four different boundary conditions mentioned above. 

We are confident that our algorithm works properly 
because the number of ground states found by our algo- 
rithm agrees with the partition function result for every 
realization and boundary condition that we have exam- 
ined. 



III. GROUND STATE DISTRIBUTION 

Before presenting the results from our algorithm, we 
first discuss the distribution of numbers of ground states 
for different realizations with varying system size L and 
antifcrromagnetic bond ratio x. All of these results 
were obtained from the partition functions of these pa- 
alizations, generated using Saul and Kardar's methodLa. 
These ground state distributions show large sample-to- 
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Boundary Condition 




a 


A 


Periodic All 
Anti-Periodic NS 
Anti-Periodic EW 
Anti-Periodic All 


3.91 ± .066 

3.94 ± .065 

3.95 ± .056 
3.90 ± .057 


1.191 ± .047 
1.199 ± .044 
1.189 ± .040 
1.201 ± .040 


589 ± 32 
588 ± 31 
591 ± 28 
591 ± 28 



TABLE I: Fits of the boundary condition distributions to a 
Gaussian with mean /x, standard deviation a, and amplitude 
A. n and a are given in terms of z, the log 10 of the number 
of ground states. Boundary condition has little effect on the 
distribution on ground states. 



sample variations for realizations with the same L and 
x. 

Figure || shows four histograms of the number of 
ground states for the four different boundary condi- 
tions of 3006 realizations of L = 10 lattices with x — 
.5. The solid linas^.an_jthe figure are fits to a log- 
normal distributionOoEJ, where the number of real- 
izations with between m and m + dm ground states is 
G(m; /i, cr)<i(log 10 (m)), with 



G(m; fx, a) 



A 



exp 



-QgioM -m) 2 

2a 2 



(2) 



Here, /i is the most probable value of log 10 (m), a de- 
scribes the width of the distribution, and A is a nor- 
malization constant. Table | gives the parameter values 
from a x 2 fit with the errors for the bin heights taken 
to be y/Nb, where Nb is the number of realizations in a 
given bin. The distributions fit a log-normal distribu- 
tion extremely well, and all the parameters of the fit for 
the four different boundary condition are consistent with 
each other within error bars. 

This log-normal distribution means that the variations 
in the number of ground states of different realizations 
are enormous. Sampling a few realizations will not give 
a meaningful picture of the behavior of the system. Av- 
erages over realizations need significant numbers to pro- 
duce reasonable results, and still may not give sufficient 
information. Also, because the distribution of ground 
states is so wide, our methods to find ground states and 
apply perturbations to them have wildly varying perfor- 
mance on realizations with the same L and x. A few 
outliers with many ground states will completely domi- 
nate the computation time of all algorithms. A change in 
thinking is necessary - the concept of an average realiza- 
tion or number of ground states is not necessary useful 
in considering the behavior of this system. 

Figure || also demonstrates that at x = 0.5 the bound- 
ary condition does not affect the ground state distribu- 
tion. All future results in this section will be presented 
for the fully periodic boundary condition. 

Next we study how the distribution of ground states 
varies with x. Figures [?] and || show how the parame- 
ters n and a characterizing the mean and the width of 
the distribution change as x is varied between and 0.5. 
Both the mean and the width tend to increase with x 
until they saturate between x = .25 and x = .3. 



The saturation of the ground state distribution at 
x = .3 appears to be completely distinct from the break- 
down of ferromagnetic order at x ~ .1, and seems to 
be relatively insensitive to changes in system size. Since 
the distribution of ground states is essentially unchanged 
from x = .3 to x = .5, this suggests that systems in this 
range of parameters have no essential physical differences. 

Finally, we present the variation of the distribution of 
ground states with L at x = 0.5. As Figure || shows, 
increasing the system size moves the ground state dis- 
tribution over to larger numbers of ground states but 
does not change the log-normal distribution of the states. 
Again we fit these distributions to the form of Eq. |[ As 
seen in Fig. the mean of the ground state distri- 
butions at x = 0.5 scales exponentially with lattice area 
L 2 ; oc e bL2 , with b 0.03. It is because b -C In 2 
that our algorithm finds all the ground states much more 
efficiently than an exhaustive search of all configurations. 

We also investigate whether the distribution of ground 
states is self averagingE^I, that is, whether the ratio of 
the width of the distribution to its mean vanishes in the 
thermodynamic limit L — > oo. Figure 11, a plot of a/fi as 
a function of L, shows that a/fi actually increases with 
L up to L — 8 (a cautionary note for higher dimensions, 
where studies of L > 8 are extremely difficult), but that 
a /n decreases as L increases for L > 8, consistent with 
self- averaging behavior for the limit L — > oo. Data for 
larger lattices would be extremely helpful in determining 
the behavior of a//i as system size increases. 

The wide distribution in numbers of ground states can 
be rationalized in terms of the matching problem by not- 
ing that random arrangements of bonds produce rela- 
tively random arrangements of frustrated plaquettes. If 
we treat the number of possible edge choices for different 
plaquettes as random variables, and assume that they are 
essentially independent from plaquette to plaquette, then 
the number of total possibilities is multiplicative, and it is 
well-known that multiplicative, random processes lead to 
log-normal distribution sr 8 H 49 H 5C I Moreover, since entropy 
is the logarithm of the number of accessible configura- 
tions, it is natural to expect log(m) to be normally dis- 
tributed. However, though these simple considerations 
make the log-normal form plausible, they do not address 
the distribution width. 



IV. GROUND STATE ALGORITHM 
PERFORMANCE 

In this section we present the performance of our 
ground state algorithm and discuss the steps limiting its 
performance. One typically investigates the performance 
of an algorithm on a finite lattice by showing its scal- 
ing behavior relative to system sizeEj. In this instance, 
such an approach does not provide much understanding 
because of the log-normal distribution of ground states 
for a given system size and x. Any average over such a 
broad distribution would not convey much information. 
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FIG. 6: Histograms of the number of ground states for 3006 realizations with L — 10 and x = .5 for the four different boundary 
conditions. The solid lines are \ 2 fits to the log-normal distribution, Eq. [| with parameters are given in Table Q. At this value 
of x, changing boundary conditions has very little effect on the distribution of ground states. 



We instead study our algorithm's performance relative 
to the number of ground states, which we believe is a 
much better measure of performance. In addition, be- 
cause the number of ground states of a realization can 
be generated quickly from the partition functiont3, this 
measure gives a useful predictor of time to completion 
before the algorithm is run. 

For all of the following results, the algorithm was run 
on a Pentium III 800 MHz machine, with 512 MB of 
RAME3. 

Figure |l2| is a plot of the algorithm's run time versus 
number of ground states. The plot also shows power law 
fits of the results for different system sizes. Though the 
scatter is substantial, it can be seen from the plot that the 



run times remain reasonably well-behaved as the number 
of ground states increases. Assuming enough memory is 
available, we expect that systems with 10 8 ground states 
could be completed in a week of running time. 

Figure |l3| demonstrates a more useful way to charac- 
terize the scaling behavior of our algorithm. Here we 
plot the run time versus the number of matchings (sets 
of edge combinations) explored. Results from different 

system sizes agree and scale as w 0(rim 5 ), where n m is 
the number of matchings. This measure of system perfor- 
mance is not a very good predictor, because it is difficult 
to determine the number of matchings before running the 
first part of the algorithm. In the fits for both Figures [l~2| 
and (0| runs with very small numbers of matches and 
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FIG. 7: Plot of the parameter jj, defined in Eq. ^| (which is the 
base- 10 logarithm of the maximum of the ground state dis- 
tribution), versus x, the fraction of antiferromagnetic bonds. 
Data for systems with L = 6,8, 10 are shown; the qualita- 
tive features do not exhibit strong L-dependence, except for 
the overall exponential scaling of the number of ground states 
with system size. 
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FIG. 8: Plot of the parameter a defined in Eq. ^, which is a 
measure of the width of the log-normal distribution of ground 
states, versus x, the fraction of antiferromagnetic bonds. Data 
for systems with L = 6,8, 10 are shown; the qualitative fea- 
tures do not exhibit strong //-dependence, except for the over- 
all exponential scaling of the number of ground states with 
system size. 



ground states are not included, because there is an ini- 
tial overhead in computation irrespective of the scale of 
the problem. 



Figures [12| and 13 show data for different lattice sizes L 
with the same x = 0.5. Figure |IJ demonstrates that the 
same scaling of algorithm run time with ground states is 
observed if we vary x and keep L fixed. 

The algorithm scales as a power law with number of 
matchings, because the core of the algorithm (its second 
part) runs sequentially through all possible matchings of 
edges in ground states. We thus expect it to scale as 




bet of ground states 



FIG. 9: Histograms of the number of ground states for real- 
izations of different system sizes L — 4, 6, 8, 10, 12, and 14 at 
x — 0.5. The solid lines are \ 2 fits 1° log normal distributions. 
The parameters /i (which describes the mean) and a (which 
describes the width) are defined in Eq. ^| they are given in 
terms the base- 10 logarithm of the number of ground states. 
Increasing system size increases both fi and a. 



0(n m ), where n m is the number of matchings. 

The algorithm's scaling with ground states is more 
complicated, precisely because the algorithm deals with 
matchings instead of ground states, except when insert- 
ing into the list of ground states, which should go like 
0(n^ c ), where n gs b c is the number of ground states for 
a given boundary condition. In practice, the number of 
ground states is smaller than the number of matchings 
by about an order of magnitude. Also, since we observe 
that ground states are split relatively evenly among the 
boundary conditions in realizations with large numbers 
of ground states, we expect n gs b c to be about one and a 
half orders of magnitude smaller than n m . Since signifi- 
cantly more work is done only on the ground state match- 
ings, which require an additional investment of searches, 
the interplay between the time for matchings and for the 
ground state conversions is nontrivial. 

The algorithm performs well only if the ground state 
matchings are a significant fraction of all possible com- 
binations of viable edges. If two distinct systems have 
the same number of ground states, but one has 10 times 
as many matchings as the other, their run time will be 
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FIG. 10: The parameter /j, of the ground state distribution, 
defined in Eq. bl versus system area L 2 for systems with x = 
0.5. The solid line on the graph is a fit to the form /i = 
(a + bL 2 ), with a = .874 ± .025 and b = .03058 ± .00023. 
The most probable number of ground states, 10 M , increases 
exponentially with L 2 . 
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FIG. 11: Plot of the ratio a/fi versus L. This is a test of 
self-averaging for the free energy. Note that a/fi increases 
initially and then begins to fall off as L increases in a roughly 
linear way. 



very different. We denote those matchings that lead to a 
ground state as ground state matchings, and the ratio of 
ground state matchings to number of total matchings of 
viable edges is inversely correlated with the run time of 
the algorithm. This is shown in Figure [l5| a plot of run 
time versus ground state matching ratio. Realizations 
with small ground state matching ratios have large run 
times, because the algorithm spends most of its time cy- 
cling through matchings that do not yield ground states. 
To illustrate the behavior of the ground state match- 



; * x = 0.5, fit slope 0.71, 300 configurations 
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a x = 0.5, fit slope 1.16, 135 configurations 
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FIG. 12: Algorithm run time versus the number of ground 
states. The fits are power law fits to each of the separate 
data sets with more than 100 ground states. 



x L=6 x=0.5, fit slope 0.74, 300 configurations 
□ L=8 x = 0.5, fit slope 0.81, 200 configurations 
a L=10 x=0.5, fit slope 0.81, 135 configurations 
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FIG. 13: Algorithm run time versus the number of match- 
ings. The power law fits yield exponents that agree and give 
a scaling of 0(n 4 ^ 5 ). Lattices with less than 100 matchings 
were excluded from the fit because of start-up time costs. 



ings, we present in Figure [l6| a histogram of the ground 
state matching ratios for different lattice sizes L. Real- 
izations with small ground state matching ratios occur 
rarely, but increase in frequency with increasing system 
size. However, even for L = 10, over 78 percent of the 
configurations have matching ratios greater than 0.1. We 
thus expect this ratio to remain appreciable for somewhat 
larger system sizes. 

The above plots do not show ground state numbers 
higher than 2 x 10 6 , because of memory limitations of our 
hardware. The code currently stores all of the ground 
states while running to prevent duplication of ground 
states. With more memory, much larger numbers of 
ground states could be investigated relatively quickly 
with this algorithm, as seen in the scaling above. In ad- 
dition, if the algorithm could be modified so that it did 
not require all the ground states to be stored, this mem- 
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- • x = 0.1, fit slope 0.54, 200 configurations 
; □ x = 0.1. fit slope 0.85, 200 configurations 

a x = 0.2, fit slope 1, 172 configurations 
^© x = 0.3, fit slope 1.04, 26 configurations 
; x=0.4, fit slope 1.04, 53 configurations 

- H x = 0.5, fit slope 1.16, 135 configurations □ 




Number of ground states 



FIG. 14: Algorithm run time as a function of the number of 
ground states for L = 10 realizations with various x values. 
The power law fits are to each of the separate data sets. 
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FIG. 15: Algorithm run time versus the fraction of all match- 
ings constructed from viable edges that are ground state 
matchings. Small ratios imply long run times. 



ory constraint would be removed. This would require 
some cleverness about how matchings are converted into 
ground states. One could conceive of keeping only a list 
of "problematic" ground states that correspond to mul- 
tiple matchings, if they could be determined easily. 



V. DESTRUCTION OF FERROMAGNETIC 
ORDER 

In this section we investigate the destruction of ferro- 
magnetic order that occurs aSpa.JJiiPiifrari.inin nf.antifer- 
romagnetic bonds, is increasedlj'HSHEpiEiEii^II. 

As Bendisch and collaborators discusst^OEll, investi- 
gating how the ground state energy depends on bound- 
ary condition is a powerful method for locating the tran- 
sition at which ferromagnetic order is destroyed. In a 
system of infinite size, when x is less than the transi- 
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0.01 



Ground state matching fraction 

FIG. 16: Histograms of the ground state matching ratios for 
realizations with L = 6, 8, and 10. In all cases, at least half 
of the realizations had ground state matching fractions equal 
to 1. Note the decreasing ground state matching fractions as 
system size increases. 



tion point x c 0.1, one expects all the lowest energy 
states to occur when the boundary conditions are con- 
sistent with ferromagnetic order, while for x > x Cl there 
should be no preference for this type of boundary con- 
dition. Refs.tEKO present calculations supporting this 
picture for different non-toroidal boundary conditions. 
Our results obtained by computing the partition func- 
tion for toroidal boundary conditions also support this 
picture. However, because of our relatively small sys- 
tem sizes, finite size effects in our calculations are large. 
We did not do a finite-size scaling analysis of our data 
because we do not expect it to yield qualitatively new in- 
formation about the destruction of ferromagnetic order. 

Our algorithm enables us to investigate in detail how 
the ground states change jwhen ferromagnetic order is 
destroyed. Barahona et al.E3 present evidence that long 
range order in the system is related to whether a set of 
spins with the same relative orientation in all the ground 
states spans the system. We can refine this picture fur- 
ther by noting that the set of ground states can natu- 
rally be subdivided into clustersEj, where a cluster is a 
group of ground states related by a sequence of single 
spin flips, each of which leaves the energy the same. By 
definition, all states in the same cluster can be reached 
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FIG. 17: Three bunches of fhppable spins. The fhppable spins 
are denoted by filled-in circles and the frustrated bonds by 
thick lines. 
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TABLE II: Table of the runs used to generate Figure |l9| N t 
is the number of distinct lattices, N r is the number of ground 
state realizations, and N scr is the number of ground state 
realizations with a single ground state cluster. 



from each other by single spin flips without raising the 
energy, whereas ground states in different clusters can 
only be reached from each other without raising the en- 
ergy by making cooperative flips of multiple spins. A 
realization's ground states may all fall into a single clus- 
ter or populate many distinct clusters. For our systems, 
the number of clusters is moderate — we have observed up 
to 12 clusters in a single 10 x 10 realization. It is natu- 
ral to ask whether the destruction of ferromagnetic order 
corresponds to growth of the number of spins contribut- 
ing to the individual clusters, or whether the relationship 
between the clusters plays a vital role. 

To address this question, it is useful to focus on the 
spatial relations between the ground states in each clus- 
ter. To do this, we first define a flippable spin as one 
with an equal number of satisfied and unsatisfied bonds; 
all the states in a cluster are related by sequential flips of 
flippable spins. We then define a bunch to be a group of 
flippable spins in a given cluster whose flippability does 
not depend on the state of other flippable spins in the 
system. Figure [I?] shows three spin bunches. 

The first bunch, spin A, is just a single isolated spin. 
The second bunch consists of spins B and C. Note that 
if spin B is flipped, then spin C is no longer flippable, 
and conversely. The third bunch consists of spins D, E, 
F, and G. The bunch contains all four spins, because if 
G is flipped, then F becomes flippable. If D and F are 
both then flipped, then E becomes flippable. 

Identification of bunches gives a complete picture of 
the clusters. Different clusters cannot have all the same 
bunches, though a given bunch can appear in more than 
one cluster. Bunches are useful because within a given 
cluster they are independent, so their contribution to the 
ground state degeneracy is multiplicative. 

We extract from the complete set of ground states all 
the bupches of a system using an algorithm described 
m Ref£3. Fi gure lq shows bunches from three realiza- 
tions with x = 0.05 < x c , x = 0.1 ~ x c , and one with 
x = 0.15 > x c . One can see that the bunch structure 
for a single cluster does not change drastically as one 



crosses the transition, but when x > x c multiple clusters 
exist and the overlap of all the bunches from the differ- 
ent clusters spans the system. We believe that the key 
element governing the destruction of ferromagnetism is 
whether the overlap between the different bunches in dif- 
ferent clusters is such that the union of all the bunches 
forms a path that percolates across the system. Thus- 
the "rigidity" transition discussed by Barahona et al.cij 
is governed by overlap of the bunches composing different 
clusters. 



Figure 19 shows the single cluster fraction as a function 
of x for L = 6, 8, and 10. The number of realizations 
used to generate this figure is shown in Table ||. The 
presence of multiple clusters is strongly correlated with 
the destruction of ferromagnetism. 

We can investigate further the bunch overlap as x is 
increased through the spin glass transition. We do this 
by defining an overlap fraction for our realizations. To do 
this, we first sum up the number of spins in bunches for 
each cluster in the realization, n s b- Then, all the bunches 
in all the clusters of a given realization are overlaid, and 
the total number of spins covered by bunches is counted, 
o s b. We define the overlap fraction Of as 



°f 



n s b — o s b 
n s b 



(3) 



If no spins overlap between the bunches in different clus- 
ters, o s b = n s b and Of = 0. If all the spins between 
different bunches in N different clusters overlapped (an 
impossibility), then Of would tend toward 1 as N — > oo. 
The overlap fraction 0/ is not defined for a realization 
with only a single cluster. 

Figure ^o| shows histograms of o/ for various x values 
for lattices with L = 8. As the number of realizations 
with multiple clusters rises, the overlap fraction for these 
clusters also increases. Not only are there more clusters, 
but the variability of the bunches in these clusters also 
increases. 

Thus, the appearance of multiple clusters for a sin- 
gle realization and the increase of overlaps between the 



11 



L = 1 x = 0.1 



L = 1 x = 0/ 




FIG. 18: Three different 10 x 10 realizations at x = 0.05, 
x — 0.1 w x c , and x = 0.5. All the bunches in each cluster 
are shown with hash marks, with different angles signifying 
different clusters. Multiple clusters allow bunches to combine 
to span the space and destroy ferromagnetic order. 



bunches of these clusters both seem to be closely cor- 
related to the destruction of ferromagnetic order. We 
believe these additional markers of the transition will aid 
in the understanding of this transition. 



-L □ L = 1 C 
x L = 8 
o L = 6 
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FIG. 19: Plot of the single cluster fraction vs x for L = 6, 8, 
and 10. Error bars are estimated as \/N r , where iV r is the 
number of realizations. At x = .05, in the spin glass phase, 
the single cluster fraction is essentially 1. As i increases and 
passes through the ferromagnetic transition, the single cluster 
fraction drops precipitously. 



VI. DISCUSSION 

This paper investigates the ground states of the two- 
dimensional ±J Edwards- Anderson spin glass. An al- 
gorithm that finds systematically all the ground states 
of this model is presented. This algorithm is used to 
investigate the nature of the destruction of long-range 
ferromagnetic order as the fraction of antifcrromagnetic 
bonds is increased. 

The running time of our algorithm for a given realiza- 
tion is found empirically to scale roughly linearly in the 
number of ground states for the range that we study (up 
to 10 7 ground states). Memory is the most serious lim- 
itation, because to avoid duplication each new ground 
state must be compared to those already found. Our 
method enables the examination of significantly larger 
systems than would be accessible by exhaustive enumer- 
ation of all states — for example, a system with L = 10 
has 2 100 ~ 10 30 configurations, but only 10 2 — 10 7 ground 
states. The advantage of our technique, is even greater 
when we introduce quantum tunnelingp], where the full 
Hamiltonian of a system with L = 10 is a 2 100 x 2 100 
matrix. 

We use our method to investigate the transition at 
which ferromagnetism is destroyed as the fraction x of 
antiferromagnetic bonds is increased. By comparing the 
"rigidity" of bonds between different ground state config- 
urations, we can obtain new insight into the percolative 
nature of the phase transition at which ferromagnetism 
disappears. 

In general, we find that many quantities exhibit large, 
non-Gaussian variability between realizations. For exam- 
ple, the number of ground states is well-described by a 
log-normal distribution, with the ratio of the width to 
the mean of the distribution increasing with system size 
L up to the relatively large value L ~ 8. This relatively 
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FIG. 20: Plots of the overlap fraction 0/ for L — 8 realiza- 
tions. As the number of multiple cluster realizations increases, 
the structure of the bunches also grows increasingly compli- 
cated as the overlap fraction increases. 



large length scale is unsettling when one keeps in mind 
the relatively small sizes thfl|t-c«tfi-be4nvestigated numer- 
ically in higher dimensionsBCJ'EJSE^I. 

The availability of complete information about all the 
ground states provides a new avenue for probing the na- 
ture of the system at low temperatures, and our work 
can be extended in many directions. We have investi- 
gated how the introduction of quantum tunneling and 
of coupling to a deformable lattice affect the low-energy 
landscape of the modeO, and the effects of other physical 
perturbations should be examined. Our work also pro- 
vides a means for stringent validation of other sampling 
methods that can be used to study larger systems. We 
also note that in the course of this work we have shown 
that one can adapt the standard matching algorithm to 
study systems with fully toroidal periodic and/or an- 
tiperiodic boundary conditions. It would be interesting 
to compare the effects of changing boundary conditions 
of this type to those in Ref.Efl. 
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APPENDIX A: DETAILED DESCRIPTION OF 
THE GROUND STATE ALGORITHM 

1. Goals of the algorithm 

These notes are a short description of our algorithm 
that generates all of the ground states of an Edwards- 
Anderson 2-D spin glass. This algorithm was designed 
with several goals in mind. 

The algorithm should be relatively simple, so that it is 
easy to understand and implement. 

The algorithm should proceed in a completely deter- 
ministic and systematic way, so there is no way to get off 
track or lose information. 

The algorithm should be exhaustive. Given infinite 
time and memory, it should be able to find all the ground 
states of an arbitrary lattice. 

The algorithm should use as much information as pos- 
sible. It should use the fact that selecting one edge au- 
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tomatically eliminates many other edges, simply because 
the two plaquettes joined by the edge cannot be linked 
by any other edges. It should also use the fact that the 
number of edges in at least one ground state is a small 
subset of the total number of edges. 

2. Short description 

The algorithm has two parts. The first is preliminary 
and the second uses the results of the first to generate all 
the ground states. 

The first part finds all the viable edges: edges that 
exist in matching solutions of the correct energy. This 
list of viable edges is then the only thing that needs to be 
studied, since the other edges do not exist in matchings 
with the correct ground state energy. 

The second part takes the list of viable edges and 
systematically goes through the possible combinations 
(those that connect each plaquette to one and only one 
other plaquette) and tests them to see if they corre- 
spond to ground states. Once all the combinations are 
exhausted, all the ground states will have been found. 

3. Detailed outline of the algorithm 

1. Given L and x, generate realization with randomly- 
placed bonds 

2. Find the partition function of this realization 

(a) gives number of ground states 

(b) gives energy of ground states 

3. Convert this lattice to a matching problem in graph 
theory 

(a) frustrated plaquettes — > nodes 

(b) path of unsatisfied bonds between frustrated 
plaquettes — > edges 

4. Construct the list of viable edges (edges that ap- 
pear in at least one ground state matching) 

(a) run blossom algorithm to find a matching so- 
lution of the correct energy. 

(b) add the edges in this matching solution to the 
list of viable edges. 

(c) build up the list of viable edges, one plaquette 
at a time, until all the edges corresponding to 
ground states are found. 

i. remove all edges in the list for the current 
plaquette from the problem. 

ii. run blossom algorithm to get a new 
matching solution. 

iii. if the new matching has the correct en- 
ergy, add the new edges to the list of vi- 
able edges, return to i. and do it again. 



iv. if it does not, there are no new edges for 
this plaquette. Go to the next plaquette. 

5. Take the list of viable edges and pair up the pla- 
quettes until all ground states are found. 

(a) sequentially go through the various combi- 
nations of edges until every ground state is 
found. 

i. pick series of edges until each plaquette is 
matched up with one and only one other 
plaquette. 

ii. if matching solution does not have the cor- 
rect energy, go to i. 

iii. convert this matching solution to a set of 
spins and bonds for each boundary condi- 
tion. If this corresponds to a ground state, 
enter it for that boundary condition only. 

iv. return to i. if some edge combinations 
remain. 

(b) number of possible combinations is manage- 
able, because each pick of an edge removes 
two plaquettes from the list that need to be 
matched. 

6. When a matching solution is found, it does not nec- 
essarily correspond to a ground state, so all ground 
states are kept in memory. 

(a) one matching can represent several ground 
states. (A simple example is a plaquette di- 
agonally separated by one from another pla- 
quette joined by an edge. This represents two 
ground states.) 

(b) multiple matchings can represent one ground 
state. (Take the example of four plaquettes in 
a + formation linked together. Many possible 
matchings represent the same ground state.) 

(c) all ground states are kept in memory and 
when a new matching solution is found, this 
is compared to ground states already found in 
the system. If it is new, it is added to the list. 

7. Once the algorithm is done, the number of ground 
states found is compared to the number found in 
step 2 as a check. 

a. Finding all the viable edges 

This part of the algorithm uses removal of edges from 
the matching problem (cutting) to find all the viable 
edges. It works by repeatedly cutting all the edges in- 
volving a specific plaquette until all matchings generated 
are higher energy. At that point, all viable matchings 
contain one of the edges already cut. 

As an example, suppose we have a 4 x 4 lattice with 
12 plaquettes and 66 edges. The arrangement of bonds 
and frustrated plaquettes is shown in Figure ElL 
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FIG. 21: 4x4 Lattice for algorithm example. Ferromagnetic 
bonds are dark lines and antiferromagnetic bonds are dashed 
lines. Frustrated plaquettes are represented by open circles 
and labelled by number. 

We run the blossom program without cutting any 
edges and get a ground state matching with edges 
(9,20,23,34,46,63). Edge 9 links plaquettes 1 and 10, 
so we put edge 9 in our list under both plaquette 1 and 
plaquette 10. We would do the same for the five other 
edges in the matching. 

We want to find all the edges that appear in ground 
state matchings that connect to plaquette 1. We thus 
cut edge 9 and feed the problem back into the blos- 
som algorithm. We get the ground state matching 
(1,23,34,49,55,63). Edge 1 links plaquettes 1 and 2, 
so we put edge 1 in our list under both plaquette 1 and 
plaquette 2. We do the same for the other five edges. 

We now cut both edge 9 and edge 1 and feed the re- 
sult into the blossom algorithm. We get out the result 
(8,20,23,34,49,56), a matching with an energy larger 
than the ground state energy. This means that every 
ground state matching contains either edge 1 or edge 9. 
We can then do the same for plaquette 2,3, etc. until we 
have all the edges that appear in ground state matchings, 
what we call viable edges. 

b. Generating all the ground states 

Now that we have the list of viable edges, we restrict 
our attention to it. All ground states can be found by 
combining these edges. We combine these edges such 
every plaquette has one and only one other plaquette 
connected to it by an edge. 

As an example, consider the 4x4 example from above. 

The set of possible edges, grouped by plaquette, looks 
like this: 

plaquette 1: (1,9) 
plaquette 2: (1,20) 
plaquette 3: (22, 23) 
plaquette 4: (22, 34) 



plaquette 5: (23,39,41,42) 
plaquette 6: (39,46,49) 
plaquette 7: (46, 52, 55) 
plaquette 8: (34,41,52,60) 
plaquette 9: (42,61,63) 
plaquette 10: (9,49,61,64) 
plaquette 11: (20,55,64,66) 
plaquette 12: (60, 63, 66) 

We want to generate a matching with six edges. Start- 
ing off is easy. We merely pick the first one in plaquette 
l's list, edge 1. Edge 1 involves plaquettes 1 and 2, so 
we skip to plaquette 3's list and select edge 22. That 
involves edges 3 and 4, so the next highest plaquette is 
5. We select edge 23, but that involves plaquette 3 and 
is forbidden, so we skip to the next one, which is 39. 
This continues until all the plaquettes are accounted for 
and we get the matching (1, 22, 39, 52, 61, 66). We would 
then convert this matching into possible ground states for 
each boundary condition. Those potential ground states 
that have consistent spin configurations and energes are 
entered into the ground state lists. 

Now we go to the next state. We remove the last two 
edges from our current state, giving us (1, 22, 39, 52). The 
second to last plaquette list we looked at was 9, so we 
choose the edge after edge 61. We get 63, which in- 
volves plaquettes 9 and 12. We then skip to plaquette 
10 and choose edge 9. That involves plaquette 1, which 
is already linked. Edge 49 involves plaquette 6 which is 
taken, and edge 61 involves plaquette 9 which is already 
taken. Only edge 64 involves the empty plaquettes 10 
and 11, so we get the new state (1, 22, 39, 52, 63, 64). 

We continue in this manner to get all of the ground 
states. If we reach the end of a list for a specific plaquette, 
we know to erase that entry in the current state and move 
to the previous plaquette list. If we reach the end of 
the list for plaquette 1, we know we have completed the 
algorithm and found all the ground states. 



4. Notes on the algorithm 

The algorithm currently is memory-limited. The blos- 
som algorithm and selection of edges for each matching 
in part 5 are both very fast. Eventually, as the number 
of ground states grows with the size of the system and 
the value of x, the available memory of the computer is 
exhausted. 

The reason every ground state is stored in memory now 
is that in certain circumstances, distinct matching solu- 
tions can correspond to the same ground state. If these 
circumstances were determined and recorded, then only 
a small subset of relevant ground states would have to 
be kept, while the others are written to a file. Whenever 
a matching solution with a "dangerous" edge or edges 
came up, it could be compared to the small subset of 
comparable matching solutions or ground states. It is 
not clear if this approach is really feasible or not. 

Finally, certain combinations of edges that automati- 
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cally lead to higher energy matching solutions constantly should greatly decrease the number of combinations of 
recur. If a list of these could be quickly generated, it viable edges that must be explored. 



* Electronic address: jwlandry@uchicago.edu; Present ad- 
dress: Sandia National Laboratory, Albuquerque, NM 
87185 

t Present address: Department of Physics, University of 
Wisconsin, 1150 University Avenue, Madison, WI 53706 

1 S. Edwards and P. Anderson, J. Phys. F 5, 965 (1975). 

2 A. Hartmann, Phys. Rev. E 59, 84 (1999). 

3 A. K. Hartmann, Europhys. Lett. 45, 746 (1999). 

4 M. Palasssini and A. Young, Phys. Rev. Lett. 83, 5126 

(1999) . 

5 F. Krzakala and O. Martin, Phys. Rev. Lett. 85, 3013 

(2000) . 

6 F. Krzakala and O. Martin, Europhys. Lett. 53, 749 

(2001) . 

7 E. Marinari, G. Parisi, F. Ricci-Tersenghi, J. Ruiz-Lorenzo, 
and F. Zuliani, J. Stat. Phys. 98 (2000). 

8 E. Marinari and G. Parisi, Phys. Rev. B 62, 11677 (2000). 
J . Lamarcq, J.-P. Bouchaud, O. Martin, and M. Mezard, 



33 
34 
35 
36 
37 



3!) 
40 
41 



c |)nd-mat/0107544 



A. K. Hartmann, J. Phys. A 33, 657 (2000). 

A. Sandvik, Europhys. Lett. 45, 745 (1999). 

A. K. Hartmann, Eur. Phys. J. B 13, 539 (2000). 

L. Saul and M. Kardar, Nuclear Physics B 432, 641 (1994). 

W. J. Cook and A. Rohe, INFORMS Journal of Computing 

11, 138 (1999). 

F. Barahona, R. Maynard, R. Rammal, and J. P. Uhry, J. 
Phys. A 15, 673 (1982). 

G. Toulouse, Commun. Phys. 2, 115 (1977). 
I. Bieche, Ph.D. thesis, Grenoble (1979). 

J. Edmonds, Canadian Jornal of Mathematics 17, 449 
(1965). 

E. Lawler, Combinatorial Optimization: Networks and 
Matroids (Holt, Rinehart, and Winston, New York, NY, 
1976). 

D. Applegate, R. Bixby, V. Chvataland, and W. Cook, 
Concordats A co de for solving traveling salesman vrob- 
lems, seet-J, UR L tittp : //www .math.princeton. edu/tsp/ 



iu M. Palassini and A. Young, J. Phys. Soc. Jpn. 69, 165 

(2000) . 

11 L. Saul and M. Kardar, Physical Review E 48, R3221 
(1993). 

12 C. de Simone, M. Diehl, M. Junger, P. Mutzel, G. Reinelt, 
and G. Rinaldi, J. Stat. Phys. 84, 1363 (1996). 

13 N. Kawashima and H. Rieger, Europhys. Lett. 39, 85 
(1997). 

14 T. Shirakura and F. Matsubara, Phys. Rev. Lett. 79, 2887 
(1997). 

15 F. Matsubara, T. Shirakura, and M. Shiomi, Phys. Rev. B 
58, R11821 (1998). 

16 A. K. Hartmann, Eur. Phys. J. B 8, 619 (1999). 

17 M. Shiomi, F. Matsubara, and T. Shirakura, J. Phys. Soc. 
Japan 69, 2798 (2000). 

18 A. Middleton, Phys. Rev. Lett. 83, 1672 (1999). 

19 A. Middleton, Phys. Rev. B 63, 060202(R) (2001), article 
no. 060202. 

20 M. Palassini and A. Young, Phys. Rev. B 60, R9919 (1999). 

21 J. Houdayer, cond-mat/0101116. 

22 C. Newman and D. Stein, Phys. Rev. Lett. 84, 3966 (2000). 

23 A. Carter, A. Bray, and M. Moore, cond-mat/0108050. 

24 A. Hartmann and A. Young, preprint cond-mat/0107308. 

25 M. Palassini and A. Young, Phys. Rev. B 63, 140408(R) 

(2001) . 

26 L. Bieche, J. P. Uhry, R. Maynard, and R. Rammal, J. 
Phys. A 13, 2553 (1980). 

27 K. Binder and A. Young, Rev Mod Phys 58, 801 (1986). 

28 N. Persky, I. Kanter, and S. Solomon, Phys. Rev. E 53, 
1212 (1996). 

29 W. Janke, Physica A 254, 164 (1998). 



N|. Hatano and J. E. Gubarnatis, Prog. Theor. Phys. Suppl. 



Concorde .html 

m M. Kac and J. Ward, Phys Rev 88, 1332 (1952). 

45 R. Potts and J. Ward, Prog. Theo. Phys. 13, 38 (1955). 

46 P. Burgoyne, J. Math. Phys. 4, 1320 (1963). 

47 M. Glasser, Am. J. Phys. 38, 1033 (1970). 

48 J. Aitchison and J. Brown, The Log-normal Distribution 
(Cambridge University Press, Cambridge, UK, 1957). 

49 E. Crow and K. Shimizu, Log-normal Distributions: The- 
ory and Application (Dekker, New York, NY, 1988). 

50 E. Limpert, W. Stahel, and M. Abbt, Bioscience 51, 341 
(2001). 

51 B. Derrida, Physica D 107, 186 (1997). 

52 S. Kirkpatrick, Phys Rev B 16, 4630 (1977). 

53 G. Grinstein, C. Jayaprakash, and M. Wortis, Phys. Rev. 
B 19, 260 (1979). 

54 J. Vannimenus, J. Maillard, and L. de Seze, J. Phys. C 12, 
4523 (1979). 

55 J. Bendisch, Physica A 245, 560 (1997). 

56 J. Bendisch and H. Trotha, Physica A 253, 428 (1998). 

57 M. Achilles, J. Bendisch, and H. Trotha, Physica A 275, 
178 (2000). 

58 A. K. Hartmann, Phys. Rev. E 63, 016106 (2001). 

59 J. W. Landry and S. N. Coppersmith (2001), manuscript 
in preparation. 

60 A. K. Hartmann, Phys Rev E 60, 5135 (1999). 

61 J. Houdayer, F. Krzakala, and O. Martin, Eur. Phys. J. B 
18, 467 (2000). 

62 H. Katzgraber, M. Palassini, and A. Young, Phys. Rev. B 
63, 184422 (2001). 

63 E. Domany, G. Hed, M. Palassini, and A. Young, cond- 
mat/0104264. 

64 J. Bogan, Honeycomb clu ster information, URL http : / / 



138, 442 (2000). 

31 K. Bhattacharya and J. Sethna, Phys. Rev. E pp. 2553- 
2562 (1998). 

32 H. Rieger, cond-mat/9705010. 



honeycomb .uchicago . edu. 



u!J The machine was,a part o 



a Beowulf Linux cluster (details 
available onlineLj). Since the code was not parallelized, 
each run was only on a single drone. 



